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Abstract 

We examine the dynamics of hard spheres and discs at high packing fractions in two and three 
dimensions, modeling the simplest systems exhibiting a glass transition. As it is well known, coop- 
erativity and dynamic heterogeneity arise as central features when approaching the glass transition 
from the liquid phase, so an understanding of their underlying physics is of great interest. Cooper- 
ativity implies a reduction of the effective degrees of freedom, and we demonstrate a simple way of 
quantification in terms of the strength and the length scale of dynamic correlations among different 
particles. These correlations are obtained for different dynamical quantities Xi(t) that are con- 
structed from single-particle displacements during some observation time t. Of particular interest 
is the dependence on t. Interestingly, for appropriately chosen Xt(t) we obtain finite cooperativity 
in the limit t — > oo. 
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The remarkable features of glass-forming liq- 
uids, so is agreed to presently, can widely be at- 
tributed to collective phenomena. We know that 
they become more and more important if we ap- 
proach the glass-transition point. It is a great 
challenge to understand collective phenomena 
because we expect them to be a kind of univer- 
sal origin of glassy behavior. 

The goal of this paper is to introduce a way 
of quantifying the degree of cooperativity by 
exploiting the correlations of two-time, single- 
particle quantities. This will be done for two- 
and three-dimensional systems thus revealing 
their very similar behavior. Naturally, these cor- 
relations have a spatial aspect, which shall be 
examined in detail. In principle, the idea is not 
new, because a considerable amount of work has 
recently been done on this subject [111 ^, ||, |j] . 
Even an experimental determination of dynam- 
ical length scales has become possible through 
multi-dimensional NMR ||| 0]. Particularly in- 
teresting, the theory of spin glasses makes pre- 
dictions about the behavior of dynamic suscepti- 
bilities and - connected to it - dynamical length 
scales, when approaching the glass transition 
Summing spatial correlations, the suscep- 
tibility shows a divergence for the analyzed spin 
model near the mode-coupling critical tempera- 
ture T+. There is some evidence, that this diver- 
gence is present in structural glasses, too [||, ^] . 

Up to now, the discussion of dynamical length 
scales has mainly been focussed on the mobility 
of particles, but not on their direction of motion. 
Various versions of spatial correlators are in use, 
e.g. (fi(0,t)fJ,(R,t)), where ^(t) represents the 
length of the total displacement during [0, t], i.e. 
the mobility, see [||]. As another example, one 
attributes to ynit) the value of one, if particle i 
is slow, and zero else [|J], see section 2. We will 
demonstrate, however, that the directional as- 
pect of motion is crucial for interparticle correla- 
tions, resulting, e.g., in a much stronger density 
dependence of dynamical length scales. 

From correlators of the above type, we can 
obtain the spatial extent as well as the overall, 
or mean, degree of cooperativity in the system's 
motion. In literature, one uses the detour via 
Fourier space and the fitting of Orstein-Zernike 
functions, to determine dynamical length scales. 
Throughout the present work, however, we will 
stay in real space, which will ease the interpre- 
tation of our data. 

We will present a treatment of the overall co- 
operativity, which, in contrast to existing work, 
makes possible a pictorial understanding of our 
results in terms of a reduction of degrees of free- 



dom. An important point is the time scale t 
defining the dynamical measurements. We hope 
to demonstrate, first, that dynamical length 
scales strongly depend on t, and, second, that 
t = T a is not a sensible choice. Furthermore we 
show that the reduction of the degrees of free- 
dom is directly related to the Haven ratio, well 
known to characterize cooperativity effects for 
the ion dynamics in ion conductors 1 1C ] . 

The organization of the paper is as follows. 
Section 1 gives the details of the performed 
simulations and introduces the main dynamical 
features via common single-particle quantities. 
Section 2 formulates our approach to quantify 
the system's overall cooperativity, which is well- 
known to be the integral of spatial correlations. 
The latter are treated in section 3, obtaining 
their strength and length scale. We conclude by 
a discussion of our results in section 4. 



1 Simulation details 

It is the advantage of a hard-sphere (HS) system 
for computer experiments, that the pair poten- 
tial 



x : Tij < R4 + Rj 
: otherwise 



forbids certain regions of configurational space, 
so that Monte Carlo steps are simply denied 
when particle overlaps occur. These computer- 
friendly 'yes/no' decisions make a Monte Carlo 
implementation of HS dynamics very efficient. 
The volume fraction 
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takes the role of temperature, which is not a 
relevant control parameter here. 

An important input parameter is the distribu- 
tion of particle sizes R it i.e. the polydispersity. 
It determines to a large extent, how amorphous 
the system is. For example, a bimodal mixture 
of spheres can be used to prevent crystallization. 
In this work, we use a continuous, gaussian dis- 
tribution of width crp and mean radius {R4) = 1, 
the latter serving as the unit of length. Particles 
of \Ri — 1 > 3crp were not used, because they 
would slow down the simulations very much. 
Former experiments show, that for 3d systems, 
crp = 10% is enough to obtain a stable amorphous 
state, i.e. lacking long-range order [11]. In the 



two-dimensional case (discs), we will work with 

op = 25%. 

Although a HS system seems rather artificial 
at first sight, there is great interest in its prop- 
erties, from both the theoretical and the exper- 
imental side. That, on the one hand, is due 
to the unbeatable simplicity, and on the other 
hand to the fact that HS are well represented 
by colloids in real life. Microscopically, colloidal 
particles perform free diffusion in their solvent, 
which is one of the reasons why we have cho- 
sen a Monte Carlo algorithm to generate the 
dynamics. Rather than integrating Newton's 
F = ma, we propagate the system according to 
the Langevin equation 

(i-i = -diV{r u ...v N ) + T7, ; , (1) 

where white noises rj^t) are directly coupled 
to the particles' positions. For very short wait- 
ing times t, there will hardly be any collision, 
so the potential term in Equ.(Q) can be ne- 
glected. The result is a free diffusion for t — ► 0, 
i.e. (r 2 (i)) w 2dD t, where d G {2,3} denotes 
the number of dimensions. All particles have 
equal masses, and their microscopic diffusion 
constants D will be the same. Although this 
kind of dynamics is convenient for simulations, 
the approach via Newton's equations (molecular 
dynamics) would lead to quite similar results. 
Naturally, the trivial short-time motion would be 
completely different from the Monte Carlo case, 
but the relevant information for longer times is 
expected to be insensitive to the microscopic dy- 
namics. Recently, this has been demonstrated 
for a Lennard- Jones type system [ |T^ ] . 

In a Monte Carlo step, we randomly choose a 
particle and try to displace it a random amount 
dx, whose distribution has the width A. Thus, 
A is the typical step length. We must take it 
as small as possible, because only in the limit 
A — > 0, we are sure to integrate Equ.(|lJ) cor- 
rectly. On the other hand, a too small A will 
reduce our simulation efficiency extremely, be- 
cause upon halving A, we need four times as 
many steps to cover the same distance. As a 
compromise, we try to achieve an acceptance 
rate of 50 percent, i.e. half of the displacements 
should result in valid moves, that is, producing 
no particle overlaps. Dependent on the pack- 
ing fraction, this yield values from A = 0.02 to 
A = 0.05, i.e. a few percent of the mean parti- 
cle radius. Comparisons to simulation runs with 
much smaller step sizes showed this choice to be 
sensible because no deviations arose, except for 
a trivial shift of time axis. 



In the simulations analyzed in this paper, we 
used relatively large systems in order to prevent 
major finite size effects. To be more specific, we 
have N =8756, 8960, 9201 and 9320 for the 
two-dimensional systems at ip 2 d —0.73, 0.75, 
0.77 and 0.78, which correspond to a box length 
of a hundred mean particle diameters. In the 
three-dimensional case, we used N =15422 and 
16307 particles filling a volume (50i?o) 3 , i.e. 25 
mean particle diameters in each direction. The 
corresponding volume fractions are tp 3 d =0.53 
and 0.56. 

To get a first impression of the system's dy- 
namics, it is most simple to calculate two- 
time, one-particle quantities. They show the 
same strong dependence on packing fraction as 
macroscopic transport quantities, like viscosity, 
when approaching the glass transition. The 
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Fig. 1: One-particle, two-time quanti- 

ties for the 3d packing fractions </? 3 d = 
50%, 53%, 56%, 57.3% and 58%, from left to 
right. The system sizes are N w 1000. (a) the 
mean squared displacement (r 2 (i)), and (b) the 
incoherent scattering function F^Q^m*, t). 

mean squared displacement 

(^(*)> s /^E( r *w- r *(°)) a \ 

and the incoherent part of the scattering func- 
tion at a given wave vector k, 

F 2 (k,t) = /l^e*(-W--(°))^ 

are the most common examples (see Fig. |]). 
The data shown in the figure above have been 
calculated in an earlier work, using small 3d sys- 
tems of N « 1000 particles. In the following, 
this data will not be used anymore. Interest- 
ingly, the one-particle quantities in Fig. [| dif- 
fer little from their counterparts in the large N 
systems, i.e. for the analyzed packing fractions 
tp-id = 0.53 and 0.56. The relaxation times r a , for 



instance, agree within 20% at <p = 0.56. Many- 
particle correlations, in contrast, have turned 
out to be very sensitive to system size, see be- 
low. Simulation runs for ip^d > 0.56 with num- 
bers of particles N > 10000 are not available at 
the moment, but the densities tp 3 d — 0.53 and 
0.56 seem to produce all the interesting features 
of a cold glass-forming liquid. 

We want to emphasize at this point that the 
one-particle quantities for the two-dimensional 
case look very similar to Fig. |], i.e. we 
find anomalous diffusion, as expressed by the 



resulting in (p W3i — 0.612 and Lp W2i — 0.815. 

In the case of a HS system, the reason for 
the reduction of mobility at high ip is quite easy 
to understand. The particles are tightly sur- 
rounded by the so called cages of next neighbors 
which to a large degree restrict their motions. 
On average, a particle feels a back-dragging 
force, which prevents its cage from being de- 



stroyed, \tiA |15p. If we define xj m) (e) =Ti(t = 



me) — r, (t = me — e) as the subsequent displace- 
ments of a tagged particle, the back-dragging 
force results in a negative value of the correla- 



slope of (r 2 (t)), and a plateau in the scattering tion (x^x. ( - m 1). With the Green-Kubo relation 
function F 2 , when going to high densities (not 



shown here). We can extract from (r 2 (t)) the 
ratio -2- of the long- and short-time diffusion 
constants. It describes the slowing down of the 
particles' long-distance motions upon increasing 
<p. The structural relaxation time r Q is calcu- 
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Fig. 2: Relaxation times r a (ip) and diffusion con- 
stants -jj- (ip) for 2d and 3d systems. The critical 
values ip c are determined by a MCT fit of the dif- 
fusion constant, i.e. -jy- ~ (tp c — ¥>)~ 7 - The results 
are (p cM = 58.7% and <^ 2d = 80.3%. 



lated according to the condition i<2(k, r Q ) = -, 
where the wave vector k = k max corresponds to 
the next-neighbor distance. The increasing r Q 
and the decreasing jy- indicate a great change of 
the dynamics when approaching a critical value 
of ip c ,3<) — 58.7% in three and tp cai = 80.3% in two 
dimensions, see Fig. ||| We obtained ip c by a fit 
of j^(<p) to a power law -§^(<p) ~ {<Pc — ¥>)~ 7 as 
predicted by mode coupling theory [ |l3[ ] . How- 
ever, the exact value of ip c should not be over- 
interpreted because fits to a Vogel-Fulcher be- 
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— lim 
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this immediately leads us to the conclusion that 
the so called cage effect is responsible for the 
slowing down of the motion, as expressed by 
D < D . 

On a longer time scale, the particles finally 
succeed in leaving their cages. Naturally, this 
can only occur if the neighbors rearrange in a 
collective way 

2 Cooperative Effects 

Figure |3| suggests the fact that a liquid near its 
glass transition possesses beautiful and highly 
non-trivial dynamics, see also [16, 17]. We 
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Fig. 3: Particle motions ri(i+50r Q )— n(t) during 
a time 50t q . The density is tp = 0.78 in a 2d- 
system. 



recognize regions of totally different behavior, 
some of which show very crowded and uncoop- 
erative motions, while others seem to act as one 



resulting in collective flows. The mobility ob- 
viously varies between different areas, which is 
commonly referred to as dynamic heterogeneity. 
How can we quantify the degree of cooperativ- 
ity in our system? It is possible to do this by 
comparing the fluctuation of a one-particle dy- 
namic quantity Xi with its many-particle equiv- 
alent ^Xj. For simplicity, let (Xi) = 0, which 
implies (J2 Xi) = 0. We can think of Xi to be the 
displacement vector Xi = v, = r,-(t + t') — r i (t) 
or its relative length X t = Svi = Vi — (vi), where 
Vi = |vj|. The direction of motion = Vjv" 1 is 
a sensible choice for Xi, too. In the case where 
interparticle correlations are lacking, the width 
of the distribution of X will be 

Var[]T^] =^Var[JQ]. 

Correlations, however, will increase Var [J2 Xi] 
while anti-correlations will do the opposite. 
From Fig. ^, it is reasonable to expect correla- 
tions, and we define 



' Vx -^Var^] 
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If our expectation is right, then 7V™ P will be 
larger than one. 

We now claim that Nx° p measures the total 
reduction of degrees of freedom caused by cor- 
relations. In the simple case of uncorrelated mo- 
tion ((XiXj) = 0), we obtain N™ e = 1, whereas 
the other extreme of totally correlated motion 
(X, = Xj) results in N?T = N. If, more gener- 
ally we have M identical variables Xi in each of 
L independent groups, N = ML, we obtain 



NZ 



= 1 



1 



N M-l 
i=l j=l 



These examples show that 7V™ op is indeed a 
reasonable quantitative measure for the degree 
of cooperativity. In real life, correlations will 
not be perfect, i.e. a hundred percent, and 
their strengths and spatial extensions will vary 
throughout the system. Hence, we should ex- 
pect N x " p to be an average or effective reduction 
factor for the degrees of freedom. 

It is important to note that N™° v = NT p {t) 
because the dynamic quantities Xi = Xi(t) are 
dependent on the time scale necessary for their 
definition. The calculation of N x x ' p turns out to 
be quite unconvenient. From one configuration, 
we only get one data point for the term £\ Xi 
in the numerator, i.e. Nx° p is not self-averaging. 
Thus, the simulation run has to be very long to 



aquire enough points for the calculation of the 
variance. A way around this obstacle is not to 
sum over the whole system, but only over lo- 
cal subsystems of n particles. This should im- 
prove statistics. In practice, at a given n, we 
randomly choose a central particle and add the 
closest (n— 1) neighbors using them as a subsys- 
tem. Repeating this procedure for a small num- 
ber of other central particles, we get some more 
subsystems of size n. 
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Fig. 4: N!7(t) at density <p 3d = 0.56, for dif- 
ferent numbers n of particles, that are summed 
over in Yli-Xi- From bottom to top: n = 
20, 100, 500, 1000, 2000, and 4000. Time is nor- 
malized tO T Q . 



Naturally a too small n will modify our results 
because we throw away some longer-ranged 
correlations, which can be important. An il- 
lustration for that is given in Fig. ^, where 
N c 5 °° ? belonging to X { = Svi = Vi - (v^) is plot- 
ted for different n at density ip 3 d = 0.56. We 
clearly see that it is necessary to take as many 
as n > 1000 particles because a major change 
can be found when decreasing n from there. In- 
terestingly, we obtained for a small system of 
N = 1066, ip 3d = 0.56 only a maximum AC P 
of nine when using n — 1066 for its calculation 
(data not shown). This is half of the value of 
N c s 7 at N = 16307, n = 1000, thus proving 
large finite size effects in many-particle corre- 
lations for the small system. Being conscious of 
this problem is especially important if one needs 
trustable numerical values for N^, e.g. for 
determining the exponent of divergence when 
cooling towards the glass transition, as is done 
in How is it possible that NZ P < 20 in the 
case tp 3 d — 0.56 although the finite-size effects 
of a too small subsystem n can be sensed even 



up to n = 1000? The reason is that particles 
are only partially correlated, which will become 
clear in the next section where we demonstrate 
the decay of correlations with interparticle dis- 
tance. Additionally, regions of fast particles are 
extended, non-compact objects (Fig. ^), so that 
we have to sum over larger subsystems to in- 
clude all their mobility correlations. 

On the other hand, however, taking n as large 
as possible is not always the best thing to do. 
We can see the reason for this at the example 
X = v taking n = N, i.e. all particles of the sys- 
tem. Then, Var \J2 v,] = because the simula- 
tion conserves the center of mass, which means 
it sets J2 v i to zero. Fortunately, our systems 
are large enough, enabling us to choose an op- 
timum value of n just between these two size 
effects. For a more thorough discussion see Fig. 
|7| and the corresponding text below. 

A word about error bars. The statistically un- 
certainty of simulation results is a consequence 
of the limited length T of the runs. If we as- 
sume the equivalence of ensemble and time av- 
erage, a quantity A can be determined up to ac- 
curacy Vary [A] ~ T _1 where the constant of 
proportionality is essentially the decay time of 
the autocorrelation function (A(t)A(0)}. We cal- 
culate Vary [A] by extrapolating its behavior for 
T' < T to T, i.e. the average over the whole 
simulation run. This is done for A = N c ™ p (t), 
only for few examples of t. The resulting errors 

SNT" = ± (Var [NT"})^ are given in the figure 
captions. 

Figure | shows A^° op = NZ° P for some densi- 
ties tp in two and three dimensions. As a func- 
tion of time, Ns°° p starts at short times with the 
value of one because the individual brownian 
motions in the microscopic regime are uncorre- 
cted. This is a trivial statement, so we do not 
have to demonstrate it for every density. For 
later times, iV™ p reaches a maximum which is 
strongly increasing with density. The following 
decay then takes some decades in time again. 
But as can be seen, a limiting value is hard to 
observe within simulation times, see below. 

Now, the behavior of 2d and 3d systems seems 
to be quite similar, although the maximum val- 
ues of A^°° p are larger in 2d. A small difference is 
the shift of the 2d maxima towards longer times. 
While in 3d they are found at approximately 
2t q , we find them in the 2d case at around 4r Q . 
This shift can be observed in other dynamical 
quantities as well. The simple reason lies in the 
different polydispersities erp =0.1 in three and 
<Tp = 0.25 in two dimensions. In 2d, the small 
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Fig. 5: [2d]:AC P 0) for n = 2000 at ip 2d = 
0.73,0.75 and 0.77 [3d]:The same for the 3d 
case, tp 3 d = 0.53 and 0.56. The statisti- 
cal error due to the finite simulation time is 
5NZ v {20T a ) = ±0.5 for tp 3d = 0.56. Again, time- 
normalization by r a . 



and - on average - fast particles cause F 2 (k max , t) 
to decay more quickly in the beginning, so we 
measure a systematically smaller r Q than in the 
3d case. If we defined i<2(k max , r Q ) = 0.01 in- 
stead of requiring i^2(k max ,r Q ) = i, this discrep- 
ancy would vanish. 

It is important to compare these results with 
N c °° p determined by the dynamical quantities 
X = v or X = n as mentioned above. For the 
interpretation of values of N c ™ p , this is essential, 
because different sensible quantities X should 
not produce totally different values of N c ™ p . Fig- 
ure |6| shows the case X — v for 2d and 3d, and 
X = n for 3d. Firstly, A^ oop and N c n °° p clearly 
display their long-time limits, which are equal 
and different from one. Secondly, N^° p still de- 
velops a maximum around t = r a , while A^ oop 
is a monotonously increasing function. Table |l| 
summarizes the main results about Nx° p - We can 
compare the maximum values of N™° p for equiv- 
alent densities (e.g. equivalent in the sense of 
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Although we find a random diffusion for every 
particle on a time scale t 3> r a as expressed by 
the diffusion law {r 2 (t)) ~ Dt, we should gener- 
ally not expect Nx° v to be one because the inter- 
particle correlations from shorter times are still 
accounted for in this quantity. This can clearly 
be seen in the following way: We decompose the 
displacement v.£ for t~^>r a into M pieces, each 
of them corresponding to a time step e—j^, i.e. 

M 




Fig. 6: [2d]: N^(t) for n = 2000 at <p 2d = 
0.73,0.75 and 0.77. [3d]: The same for the 3d 
case, ip 3d = 0.53 and 0.56, n = 2000 (+). Ad- 
ditionally, we see N^(t) (O). The statistical er- 
ror is <«V™ op (10t q ) = ±0.5 at ip 3d = 0.56, and 
SN c v oop (50T a ) = ±1.3 at (p 2 d = 0.77, for example. 
Time is normalized by r Q . 
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equal -tj-) and find a somewhat higher cooper- 
ativity in 2d. The absolute values of N™° ? for 
the different X agree to a reasonable extent, so 
that we are indeed allowed to interpret them in 
the sense of a reduction factor for the degrees of 
freedom. 

Let us now turn to the limiting value N™^. 



for particle i. For simplicity, let us say that inter- 
particle correlations are neglible for different 



time intervals, i.e. (A£ 
In this case, we obtain 



A™ } = if m ^ 



where we exploited time translational invari- 
ance, i.e. having an equilibrium liquid. This 
quantity, however, does not depend on time 
t = Me anymore, if e is fixed. We thus get an 
idea how it is possible that correlations persist 
for t — > oo. 

Finally, other choices of X, are possible, e.g. 
more exotic quantities like 



X, 



slow 
fast 

otherwise 



where the exact definition of fast and slow is of 
no importance as long as it is done in a sensible 
way. Such an analysis has been presented in ij^J] 
for a Lennard- Jones fluid using a 'dynamic sus- 
ceptibility' x ss instead of Nx° p - Its definition is 
quite similar to A r ™ p , measuring fluctuations of 
a many-particle, 'macroscopic' dynamic quantity 
Qss = J2 w * : 



A'S: 



N 2 



Unlike A r ™ p , a quantitative interpretation of the 
value of Xss is not obvious. 

Let us return to the role of the size n of the 
subsystems that are used for the calculation of 
Nx° p . As we said above, the subsystem should be 
large enough to include most of the long-ranged 
correlations of its particles, i.e. to reduce sur- 
face effects. On the other hand, n = N leads to 
lycoop _ q Decause f the center-of-mass correc- 
tion (cmc) . It is evident that our results will be 
influenced by the cmc even if we use n < N, 
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Fig. 7: The role of the subsystem size n in the 
calculation of N™ ^ (+), for the systems ip 2 d = 
77% (N = 9201) and ip 3 d = 55% (N = 16307). 
The solid lines are fits of the form A^(l - §), 
where the parameter N™°? n = 35 and 16.5, in the 
2d and 3d case, respectively. The □ marks the 
choice n = 2000. 



say n — 0.9N. Despite this fact, we need a 
cmc because the motion of the whole simulation 
box leads to unphysical results for N c ™ p . Inter- 
estingly, the center of mass performs a random 
walk with speed (v c 2 m (i)) = j^2dD t, indepen- 
dent of the packing fraction, which is the conse- 
quence of actio = reactio in a stochastic sense. 
Now, the subtle point is that the cm motion con- 
sists of two contribitions, firstly the random dis- 
placement U of the simulation box as a whole 
which would vanish if we embedded the simu- 
lation box in macroscopic system. Secondly, the 
random rearrangements w, of particles in our 
system that produce a contribution to the cm 
motion even if we forbid an overall drift of the 
box. Obviously, we should keep the second and 
discard the first contribution because the latter 
is an artefact of the limited system size. A cmc, 
however, will remove both. In the following, we 
estimate the resulting error in the calculation of 
jycoop -p^g uncorr ected displacement of particle 
i is Vj = w, + U, so that the cm motion becomes 
v cm = A J2k=i w fe + U, where the first term does 
generally not vanish. Calculating the numerator 
of N c v °° p in Equ.(§), we obtain 
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where, in the final step, correlations between 
the ^ different subsystems have been neglected. 
Thus', 
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is the result, which is too small by a factor of 
(1 — -S>). As is demonstrated in Fig. fj for the 
value of A^ oop (oo), this behavior can indeed be 
observed in our simulations. The renormalized 
value of cooperativity N™°? en , comes out as a fit 
parameter to the form of Equ.(||), see table [l] in 
the brackets. Additionally, Fig. fj tells us what 
choice of n is advisable because has to be in 
a region where the above fit works. If this is 
the case, the subsystem must have only small 
remaining correlations with the other N~n par- 
ticles because increasing n only reduces A"™ op 
by the 'trivial' factor of (1 — jt). The choice 
n = 2000 is marked by □ in Fig. [7]. 

In our case, at ip 2 d — 0.77, we observe a max- 
imum reduction of the total degrees of freedom 
by a factor of A^°° P n = 35. This, however, is yet 
only a moderately high density (<^ c2d « 0.8), so 
we should expect large collective effects at the 
glass transition. 

Finally we mention an interesting relation be- 
tween N™^ and the Haven ratio, relating the 
ratio of the self-diffusion constant and the con- 
ductivity in ionically conducting materials. Its 
zero-frequency limit H(0) is given by 

( ' *X^/o"*<v<(Q)vj(t)>- ^ 

Since for large times (v?(t)) = 2dDt = 
It J °° dt(v,(0)vj(t)) and correspondingly 
(vi(t)vj-(t)) = 2t J °° dt{vi(0)vj(t)}, it is obvious 
that 

kz = #(o) -1 . (5) 

Hence, as a side product we have obtained a 
quantitative interpretation of the inverse Haven 
ratio as the reduction of the effective degrees of 
freedom. 



3 Spatial Correlations 

The snapshot of the dynamics in 2d (Fig. |||) 
demonstrates that large spatial correlations are 



present in our systems. In the following, we 
want to quantify them as a function of the 
timescale of dynamics. 

As in the treatment of N™ *, a dynamical vari- 
able Xi should be given for each particle, again 
with the restriction (Xi) = 0. A spatial correla- 
tor can then be defined by 

(X(O)X(R)) = /^Y^XiXjSin-in-r^ . 

(6) 

Again, Xi denotes a dynamical quantity con- 
nected to the motion of particle i during the time 
interval [t a ,t +t]. Because of symmetry reasons, 
the positions r, = rj(to + f ) are used. Averaging 
over the solid angle of R, i.e. 

1 



(X(0)X(R)) 



AnR 2 



dQ(X(0)X(R)), 



results in a loss of information because the 
direction of motion of particle i, for example, 
breaks the isotropy. Figure |§] explains this in 
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Fig. 8: Correlations of the direction of displace- 
ments, Xi = n, at ip = 77%, t = I0r a . The large 
arrow in the middle shows the direction of mo- 
tion of the reference particle. For details, see 
text. 



a pictorial way: particles "in front of" or "be- 
hind" the central one have a very long-ranged 
directional correlation, while perpendicular to 
the direction of motion, we observe a kind of 
back-flow behavior which is well-known from 
[jlS|]. How has this plot been produced? First, 
we calculate all Xi — n,, i.e. the directions of 
displacements during a time interval [t ;t + t]. 



We then choose particle i and turn the whole 
(2d) system so that X[ points along the posi- 
tive x-axis. Now, the directions Xj are added to 
the average at the positions r'j (to + 1). As result, 
we obtain the field (X'(r')), which for large | |r'| | 
consists of very short vectors. Hence, for reasons 
of visualization, we plot the normalized version 
of (X'(r')) in Fig. |. 

Being aware of the complicated behavior in 
Fig. ||, let us for the moment and for simplicity 
ignore the angle-dependence, and treat correla- 
tions only as a function of interparticle distance 
R. We define the dimensionless quantity 



_ (X(0)X(R)) 
Sx(R,t)= (X2(Q)) , 



(7) 



where again its dependence on time-scale t 
should be kept in mind, just as in the case of 
NT 9 - Possible choices are Xi = Svt, Xi = v, 
or Xi = rii, where = v.i(t) - r;(0), Sv z = 
Vi — (vi) and = Vj. The functions Ss v (R,t) 
and S n (R, t) count correlations of both slow and 
fast particles because both sorts are weighted 
similarly. To be more precise, the slow particles 
are not suppressed as in the case of S V (R, t). 

Calculating Sx (R, t), we encounter a problem 
that is related to the system size. If we want the 
system's center of mass to be constant, we have 
to correct the particles' motions. But this intro- 
duces an anti-correlation of two formerly uncor- 
rected particles. As a consequence, Sx(R,t) 
will approach a negative value for large R, in- 
stead of zero. Without center-of-mass correc- 
tion, we would measure a positive number in 
this limit because the center of mass performs 
a slow diffusive motion, see the previous sec- 
tion. For the present systems, these offsets were 
smaller than 0.002, which is a small fraction of 
the amplitude of correlation. Hence, their sub- 
traction from Sx(R,t) left the function nearly 
unchanged. 

Figures || and |To| show S Sv (R,t) for the 2d 
systems f 2 d = 73% and ip 2 d — 0.77. Obviously, 
Ssv(R, t) can be described by an exponential 



S Sv (R,t) w A(t)exp 



R 



(8) 



to a good approximation, if R > 5. It is impor- 
tant to note that the amplitude A(t) is not neces- 
sarily equal to one as suggested by the definition 
of Sx(R,t) for R — > 0. In other words, the ex- 
trapolation of complicated inter-particle corre- 
lations to the one-particle quantity Sx(R = 0, t) 
would be unphysical. 




We find large deviations from the exponential 
at distances R < 5. This can be understood qual- 
itatively because certain information about the 
local packing is available. For instance, R < 1 
can only occur for two very small particles (re- 
member (Ri) = 1) which on average are much 
faster than the others. Thus, (Sv(0)Sv(R)) will 
be quite large. The oscillations for R < 6 must 
have a similar reason, i.e. special local pack- 
ings that are favorable or not for the value of 
(Sv(0)6v(R)). We can imagine that for larger R, 
the possibilities of packing become so many that 
they average Sx (R) to a structureless exponen- 
tial. This is the case for the structure factor g(R), 
too. 

Figure 11 shows Ss v (R,t) for the three- 
dimensional case. The situation is quite the 
same as in 2d, except for the long-range os- 
cillations, especially of Sg v (R : t — r a ). They, 
too, indicate the existence of structures that are 
favourable for dynamical correlations. For 2d 
systems of smaller polydispersity ap = 10%, 
which are not shown here, we find the same 
oscillations. In this case, they could be proved 
to result from local crystalline order, which oc- 
curs, if - by coincedence - many particles of ap- 
proximately the same size come together. Al- 



though the system is in an overall amorphous 
state, the small polydispersity makes local crys- 
tallites more probable, thus creating regions of 
low mobility. The oscillations for three dimen- 
sions, however, are not understood yet. 

In any event, we can extract from Ss v (R,t) 
the amplitude A(t) and the correlation length 
£x(t) as a function of the dynamic time scale 
t. The simulation runs, by the way, have to 
be much longer than the maximum time scale 
shown because the functions S$ v (R, t) are quite 
demanding with respect to statistics. For in- 
stance, distant particles, which are uncorre- 
cted, have to average (X(O)X(R)} to zero. The 
statistics M enters by a factor of -4=, so an im- 
provement of the result has a high price. Ad- 
ditionally, the dynamic heterogeneities, as visu- 
alized in Fig. ||, are very long-lived JT^], i.e. 
possess typical life-times of tens to hundreds of 
T a , dependent on tp. Thus, if we want to aver- 
age over different dynamical situations, we need 
data for some hundreds of r Q . Table ||] shows the 
lengths of our simulation runs in units of r a for 
the analyzed 2d and 3d packing fractions, re- 
spectively, including the system size N. 

Let us turn to another dynamical quantity 
Xi — rii, the direction of displacement. As 
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we see in Fig. 12, S n (R,t) is quite similar to 
Ss v (R,t), i.e. we find an overall exponential 
decay of correlations with distance R. Its char- 
acteristic length i n (t), however, is much larger 
than the previous £ 4 „(t) (see below). Again, the 
situation is quite the same in three dimensions 
(Fig. 13), i.e. exeeds £ 5 „ for ip 3d = 0.56. A 
remarkable difference to the 2d case are the ex- 
treme oscillations of S n (R,t) for t — j^r a . This 
is not understood yet, but could be due to the 
lower polydispersity in the three-dimensional 
system causing locally less amorphous packings. 
The vectorial correlation S v (R, t) is also shown 
in Fig. O, where it seems that the amplitude 
A v (t) is a constant for all time scales t > r a . 



Figure 14 summarizes the data for t;x(t), 
Xi = Svi. First of all, we notice an increase 
of correlation lengths with density. For a 3d 
Lennard -Jones system [|lp and polymers [0], 
this has already been demonstrated for the spe- 
cial choice of time scale t w r Q . But for longer 
times, even larger correlation lengths can be ob- 
served as shown in the figure. Interestingly, 
(,x{t) is a monotonously increasing function, 
with a limiting value £x (oo) = lim t _>oo £x (t) ■ For 
comparison, Fig. 14 shows the length scales for 
the directional correlation £ n (£) at ip 2 d = 0.77 
and tp 3d = 0.56. 

Table ^ summarizes the data for £x(oo), 
where the error due to fitting the exponential 
is less than ten percent. (The quantity R €UIl will 
be explained below.) The statistical error due 
to finite time averages is small enough to be 
included therein. As we see, the length scale 
£s v (po) takes a snuggish rise, growing from 4.1 to 
7.15 between ip 2 d = 0.73 and 0.77. In contrast, 
£ n and £ v (oo), starting from about the same ini- 
tial value, end up at a value of more than twice 
£^(00). The underlying physics of this very dif- 
ferent behavior of mobility and directional cor- 
relations is unclear. As in the case of Nx° p , the 
vectorial quantities rii and v^ show a very simi- 
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Fig. 13: Spatial correlation S n (R; t) and S V (R; t) 
at ipzd = 56%. 



lar behavior of their inter-particle correlators. 

So we can say, that the overall cooperativity is 
determined by both its length scale £x(t) and its 
strength, or amplitude, A(t). Eqns. @ and (0) 
show that 
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dRp(R)S x (R), 



where p(R) is the average number of particles 
found at distance R from a particle at the ori- 
gin. In a way, this is a trivial result because the 
sum over all spatial correlations should be the 
overall cooperativity. When the approximation 
Sx{R) ~ Ae~ R ^ is valid and at homogeneous 
density, we find NT"(t) ~ £ x (t) 3 A(t). The de- 
viations from the exponential for small R can 
modify this argumentation, but are unlikely to 
totally change the picture. 

Consequently, the results for iVa° op (i) and £g v 
are only compatible, if the strength of correla- 
tion As v (t) will tend to zero for long times. We 
can observe the decrease of Ag v (t) clearly in Fig. 
[9| 0/J = 0.73), but it is harder to see at higher 
densities because of the limited time window. 
For Xi — Vj, in contrast, we need a limiting 
value of A v greater than zero, if S v (R, t) is to 
be compatible with N^°° p (t) for t — > 00. Figure 
13 proves this to be the case because A v (t) is 
constant for t > r a . 

Let us finally return to the detailed picture 
of Fig. |8| It suggests, that spatial correlations 
along the direction of motion will be very differ- 
ent from them perpendicular to it. We can test 
this by restricting the summation in Equ.(^|) to 
certain angles ijj between the motion of particle 
i and the connection vector Vj— Vi. For example, 
the condition ip e [f , \ ] chooses only particles 
j that are collateral to particle i with respect to 
its motion Vj. To select in direction of motion, 
we demand ip e [0, ^] U [^§7r, it], for example. 
The restricted sum is then carried out to obtain 
S^(R, t) or SH(R, t), see Fig. |l|. Because of the 
backflow, we expect sj](i?, i) to become nega- 
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dRp{R)(X(0)X(R)) 



tive for large R > 20. This can be observed in 
the figure. In contrast, Sn(i?, t) approaches zero 
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Fig. 14: [2d]: The dynamical length scales 
£s v (t) at packing fractions ip 2 d = 73%, 75%, 77% 
and 78%, from botton to top (+). For compar- 
ison: £ n (t) (□) and ii? curl (t) (O) at ip 2d = 77%. 
[3d]: Again £tv(t) f° r the 3d densities ip 3 d — 
53% and 56% (+) and £ n (f) at <p 3d = 56% (□). 
Errors due to fitting are estimated to be less than 
ten percent. 



f o R — > oo. We notice the difference of length 
scales and £„, which is summarized in Ta- 
ble ||. The estimated error of these length scales 

Table 4: 



10 15 20 25 30 35 40 



Fig. 15: Spatial correlation S n (R,t) at tp 2 d = 
77%, t= 10r a and <^ 3(i = 56%, t = r a , calculated 
dependent on the angle ip = (v<, r^— r^). Paral- 
lel (||) means ?/> S [0, ^] U [^7r,7r], perpendicu- 
lar (_L) means tp e [|, ^p], and bulk stands for 
V; € [0,7r]. In the perpendicular case, S n (R,t) 
becomes zero at R rs i? curl and is negative for 
larger distances, so that it is necessary to plot its 
absolute value. 
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is less than 20 percent, but the exact values are 
not of interest here. Instead, the notion of a very 
different behavior of particles parallel and per- 
pendicular to the motion of the central one is 
justified for three as well as for two-dimensional 



The time dependence £,x{t), if at first sight 
surprising, can be understood quite pictorially. 
Without any further information, the probability 
for the motion of some tagged particle i is equal 
in every direction. However, if we know that, in 
the meanwhile, one or more of its next neigh- 
bors perform some specified displacements, this 
will influence the probability of movement of 
the tagged particle. More distant neighbors will 
do this as well, but the information about their 
motions has to be 'submitted' to particle i via 
nearer neighbors. It is not hard to imagine 
that the information spread can only take place 
with a finite velocity. Thus, short-time motions 
will be accompanied by a reaction of few neigh- 
bors, while long-lasting displacements will in- 
volve many of them. The monotonously grow- 
ing length scale of dynamic correlations is the 
natural consequence. 



In the limit t — > oo, we can argue as in the case 
of Nx° p - For long waiting times, the displace- 
ments u 4 = r<(|) - r;(0) and w< = r,(t) ~ r,-(|) 
become independent to a good approximation. 
Thus, correlations on time scale t can be ex- 
pressed through Ui and w t , using v, = + wj: 

(v(O)v(i?)) « (u(0)u( J R))' + (w(O)w(i?))" . 

This results in 

S v (i?,t)^<(i?,0+i<(i?,£). 

The primed and double primed versions of 
S V (R, t), respectively, denote measuring the dis- 
tances R at the end or at the beginning of the 
time interval [0, t]. On the other hand, S v (R,t) 
is defined by using the interparticle distance R 
in the middle of this interval. However, the defi- 
nitions S V (R, t), S V (R, t) and S v (R, t) produced 
the same results in our simulations, which is 
not shown here. Thus, in the long-time limit, 
the spatial correlations of the vectorial displace- 
ments Xi = v s ; become time-independent, i.e. 

S v (R,t) ^ S v (r,^\. 

The phenomenon of a growing dynamical 
length scale is exibited by much simpler systems, 
like a one-dimensional (closed) chain of N diffu- 
sive particles which are connected by harmonic 
springs, as described by the Langevin equation 

where the r\ n are independent white noises. Let 
us assume that N is a large number, say N > 
1000. The analytic solution of this many-particle 
problem is possible with the help of discrete 
Fourier transform. This enables us to calculate 
the displacement-displacement correlation, but 
this is not shown here. In this simple model, the 
length of correlation increases with time, too. 
In contrast to our simulations, £ v (*) grows un- 
til it has reached the system size. Only from 
this time on center-of-mass diffusion prevails, 
for which correlations among different particles 
are no longer relevant. Stated differently, apart 
from finite size effects the chain model possesses 
£ v (oo) = oo. In turn, we see the reason for a fi- 
nite value of i n our HS systems: 

"Particles simply can go out of each other's 
way". 

In other words, a particle that travels a long 
distance does not have to pull the whole system 



with it because rearrangements are possible. On 
average, this results in the back-flow behavior of 
Fig. |8| We are therefore tempted to relate the 
length scales (oo) to an inherent length of the 
back-flow pattern for long times. From Fig. ^, 
we see that the distance R CUIl (t) from the vortices 
to the central particle is the only sensible choice. 
Tabled shows the limiting values _R curl (oo) for the 
2d systems under investigation. Interestingly, 
i? curi (oo) is twice the correlation length £ v (oo) or 
£ n . For £s v (pd), no such relation seems to exist. 

4 Discussion 

We presented detailed information about dis- 
placement correlations, which turned out to be 
of the same nature for two-dimensional discs 
and three-dimensional hard spheres. 

Using the quantities Xj = 5v iy rii and v, as 
input for NT p (t) and S X (R, t), we were able to 
measure the total reduction of degrees of free- 
dom and the spatial extent of correlations, re- 
spectively. The data of 7V™ P is found to agree 
for these three choices of Xi, supporting the 
notion of a reduced dimensionality of motion 
in high-dimensional configuration space. The 
length scale of correlations, however, turns out 
to increase much faster with density for Xi = 
and v, than for Xi — 5vi. An explanation for this 
is lacking at the moment. Finally, we demon- 
strated, that for 2d as well as for 3d systems, an 
angle-resolved calculation of correlations is ap- 
propriate, yielding much larger length scales in 
the direction of motion than perpendicular to it. 

The important question arises how the dy- 
namical length scales £x(oo) are connected to 
static correlations. Naive attempts, using, e.g. 
spatial density correlations, have not revealed 
any significant connection. What makes things 
more complicated is £x{t)'s dependence on the 
definition of X t . On the other hand, 7V™° p (i) 
and Sx{R;t) are strongly averaged quantities, 
obtained by including many different dynamical 
situations. Thus, we should not hope to get too 
specific information from them. A deeper under- 
standing of cooperative effects will only become 
possible by a more detailed, less averaged treat- 
ment. In any way, the relation of structure to dy- 
namics is the central problem to be solved. The 
present work may help to formulate the relevant 
questions somewhat clearer. 
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